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Abstract 

The dynamics of strongly immunizing childhood infections is still not well understood. Although 
reports of successful modeling of several incidence data records can be found in the literature, the 
key determinants of the observed temporal patterns have not been clearly identified. In particu- 
lar, different models of immunity waning and degree of protection applied to disease and vaccine 
induced immunity have been debated in the literature on pertussis. Here we study the effect of 
disease acquired immunity on the long term patterns of pertussis prevalence. We compare five min- 
imal models, all of which are stochastic, seasonally forced, well-mixed models of infection based on 
susceptible-infective-recovered dynamics in a closed population. These models reflect different as- 
sumptions about the immune response of naive hosts, namely total permanent immunity, immunity 
waning, immunity waning together with immunity boosting, reinfection of recovered, and repeat 
infection after partial immunity waning. The power spectra of the output prevalence time series 
characterize the long term dynamics of the models. For epidemiological parameters consistent with 
published data for pertussis, the power spectra show quantitative and even qualitative differences 
that can be used to test their assumptions by comparison with ensembles of several decades long 
pre-vaccination data records. We illustrate this strategy on two publicly available historical data 
sets. 

Keywords: pertussis, childhood diseases, recurrent epidemics, stochastic fluctuations, power spec- 
tra 
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I. INTRODUCTION 



Childhood infections remain a pubhc health as well as a modeling challenge, despite 
many decades of control measures and theoretical efforts. Large vaccination programmes 
against some of these diseases started in the 1940s-1960s in the developed countries and led 
to dramatic reductions of incidence levels, but in the developing countries they are still a 
cause of significant levels of infant mortality IliJei. The resurgence of pertussis, also known 
as whooping cough, reported in developed countries after decades of high vaccine coverage 
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as well as a recent upsurge of measles in eastern and southern Africa prompted a 



renewed interest in assessing the efficacy of control policies for these childhood infections. 

Such an assessment must rely on simulations based on mathematical models. The com- 
plexity and diversity of the long term dynamics of childhood diseases has been long ac- 
knowledged as a major problem in mathematical epidemiology 
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15| . More recent work 



focused on contact structure, stochasticity and seasonality has brought considerable ad- 
vances in understanding and selecting some of the fundamental ingredients that drive the 
observed incidence temporal patterns 1614231] and in developing analytic tools to deal with 



these ingredients and their interplay with the model's nonlinear it ies 



2J-|32|. 



Models of higher complexity involve many parameters, some of which are difficult to 



determine, leaving considerable room for data fitting 



33|-|36|. On the other hand, one of the 



features of childhood infections is the variability exhibited by different data records, even 
within those that correspond to a single disease in comparable social environments js], isljls^]- 
Therefore, successful modeling of particular sets of incidence data records reported in the 



literature 
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39| has not closed the discussion about the key determinants of the observed 



dynamics. Measles is an exception for which a parsimonious model has been shown to 



reproduce the behavior of different data sets 
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4l[ |. and it is generally accepted that 



the disease dynamics in large urban populations is adequately described by a stochastic 
seasonally forced well-mixed susceptible-infective-recovered (SIR) or susceptible-exposed- 
infective-recovered based model. 

At the opposite extreme, pertussis keeps defying mathematical modeling, as illustrated 
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by recent efforts in different directions ll|, |39|, |42N45| . In particular, several proposals 



11 



45|, relying 



explore hypothesis about disease induced and vaccine induced immunity 
on assumptions about vaccine uptake and efficacy for the purpose of comparison with real 
data. This is a difficult problem of great interest and enormous relevance for public health. 
However, we are still lacking a sound uncontroverted model for pre-vaccination pertussis. 

Here we have sought to contribute to the goal of using pre-vaccination data records to 
obtain information about the properties of disease induced immunity. The paper focuses on 
exploring the influence of naive hosts immune response in the long term patterns of pertussis 
prevalence as given by the averaged power spectra of simulated time series corresponding 
to several decades. The power spectrum of the stochastic seasonally forced well-mixed 
SIR model is compared with those of four modifications of the model that reflect different 
assumptions about disease induced immunity. These four different assumptions have been 
proposed in the literature in the framework of deterministic models that have all been shown 
to be compatible with available data. The five variants we consider deliberately avoid all the 
complications related with contact structure and spatial spread, as well as vaccine coverage 
and efficacy and waning of vaccine-induced immunity, in order to keep a relatively low 
number of free parameters in the model. 

We show that the stochastic versions of the SIR model and its four variants have sig- 
nificantly different properties, which translate into quantitatively different prevalence and 
incidence power spectra. This opens the possibility of using the stochastic properties of long, 
well resolved data records to constrain these and other variants of the model, using the power 
spectra of the pertussis incidence time series in large urban centers in the pre-vaccination 
era as the target long term dynamics that the model should reproduce. We illustrate this 
strategy by applying it to two publicly available historical data records for pre-vaccination 
pertussis incidence. 
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II. MODELS AND METHODS 



A. Models 

We consider five seasonally forced stochastic compartmental models summarized in Fig. [1] 
as candidates for the description of pre- vaccination pertussis (the deterministic counterparts 
of these models are described in the Supplementary Material). In all cases, the population 
includes three classes of individuals, the susceptible [5* in (a)-(d), 5*1 and 5*2 in (e)], the 
infectious [/ in (a)-(c), Ji and I2 in (d)-(e)] and the recovered [R in (a)-(b) and (d)-(e), R 
and W in (c)]. Also in all five cases, there is replenishment through births of susceptible 
individuals that have never before contracted the disease. The birth rate is kept constant 
and equals the death rate yU.. 

1. SIR model 

Model (a) is the standard SIR model: recovery confers permanent immunity, and all 
infections are first infections. Infected individuals recover at a constant rate 5, and suscep- 
tible individuals are infected at a rate \t = P{t)I/N, where is the population size, / is 
the number of infected individuals and /3(t) is a periodic function of period one year that 
represents the variable contact rate associated with the school terms. For the purpose of 
studying the long term temporal patterns of the fluctuations of this stochastic process the 
particular form of is not crucial, and we shall take (3(t) = /3o(l + /3i cos27rt). 

2. SIRS model 

Model (b) is the standard SIRS model. It differs from the SIR model in that recovery 
does not confer permanent immunity. Instead, recovered individuals' immunity wanes at a 
constant rate 7 and then they become susceptible again. The fraction of primary infections 
is given approximately by fi/['~f{l — s* — i*)+fi\, where s* and i* are the equilibrium values for 
the susceptibles and infectives in the deterministic counterpart of the SIRS model (see the 
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(a) SIR 



(b) SIRS 




(c) SIRWS 



(d) SIRIS 



(e) SIRSI 

FIG. 1: The diagrams of five compartmental models for pertussis. The models (a) and (b) are 
the seasonally forced SIR and SIRS models, respectively, while the models (c), (d) and (e) are 
extensions of the SIR(S) paradigm. Model (c) allows for boosting of immunity proportional to, and 
potentially greater than, the force of infection, while model (d) allows for reinfection of recovered 
individuals. Model (e) accounts for a loss of infection-derived immunity and subsequent reinfection. 
In all diagrams, the transitions corresponding to the birth and death of individuals are not shown 
for simplicity. The meaning of the symbols associated with each transition is defined in the text. 

Supplementary Material for details on the deterministic SIRS model and on this calculation). 
Although there is no direct correspondence between primary infections and age groups in 
the model, it is reasonable to admit that a fraction of the infectives close to the fraction of 
non-primary infections, given by 7(1 — s* — i*) /[-/{l — s* — i*) + fi], are not school children, and 
that Pi should be reduced proportionally in this version of the model. This correction could 
be refined, but as we shall see in Section IIIII the properties of the stochastic fluctuations in 
this model are largely insensitive to the value of 
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3. SIRWS model 



Model (c) is an SIRS-type model allowing for immune boosting of recovered individuals 
upon reexposure to infectious individuals. In this model, infectives I recover at rate 5 and 
susceptibles S get infected at rate Aj as in the SIR model. However, in addition to the S and 
I classes, there are two classes of recovered individuals denoted by R and W . The immunity 
of the former is not permanent and they move to the waning class at a constant rate 27. 
The individuals in class W undergo two possible transitions: further immunity loss at rate 
27 to become susceptible S", or immunity boosting upon contact with infectious individuals 
to return to the recovered class i? at a rate kAc = nP^I /N , where k is the immunity boosting 
coefficient. We call this scheme the immune boosting model and denote it by SIRWS. 

An age-structured version of this model with vaccination has been used to explain the 
recent reemergence of pertussis cases despite high vaccine coverage in Massachusetts, and 



also the shifts in total and age-specific incidences before and after mass vaccination 



4-- SIRIS model 



Model (d) proposed in j45| sets a scenario based on the SIRS model with a moderate rate 
7 of immunity loss where recovered individuals are immune to severe disease but susceptible 
with reduced susceptibility to mild forms of the disease. The classes of individuals infected 
with severe and mild infections are denoted by Ji and J2, respectively. Recovery from both 
forms of the disease takes place at the same rate 5. Infectiousness of mild infections is 
reduced by a factor 77 e [0, 1] and susceptibility of recovered individuals to mild infections 
is also reduced by a factor a G [0, 1] with respect to the susceptible. Moreover, since mild 
infections typically occur in adults that are not affected by seasonal forcing the force of 
infection is therefore taken to be \t = {f3(t)Ii + f3oril2)/N for susceptible individuals S and 
Ac = (^o{h + V^2)/N for recovered individuals R. We call this scheme the reinfection model 
and denote it by SIRIS. 



The main feature of this model is that it exhibits a reinfection threshold 



46| . that is, a 
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value of infectiousness above which total disease incidence rises by one or more orders of 
magnitude, due to high incidence levels of mild infections. These mild infections are sub- 
clinical and so the disease incidence/prevalence and density of infectives in equilibrium for 
the deterministic counterpart of the SIRIS model is associated with the class Ii. 



5. SIRSI model 



Finally, model (e) proposed in 42| assumes an immune response that is a combination 
of (b) and (d), in the sense that recovered individuals are fully immune to the disease, 
but they lose this immunity at a certain rate 7 to become susceptible to repeat infections, 
although with reduced susceptibility. The two classes of susceptible individuals are denoted 
by 5*1 for the naive susceptibles and 5*2 for the susceptibles generated by immunity waning, 
and cr G [0, 1] is the reduced susceptibility factor for repeat infections. The classes of 
individuals infected with first and repeat infections are denoted by Ji and I2, respectively. 
Recovery from both classes takes place at the same rate 6. The class of repeatedly infected 
individuals contributes with reduced infectiousness to the pool of infectives responsible for 
disease transmission, and r/ G [0, 1] is the reduced infectiousness factor. Repeat infections 
play a similar role in this model to that of mild infections in the SIRIS model, and the force 
of infection is taken to be = + f3oril2)/N for Si and Ac = f^oih + V^2)/N for 5*2, 

because repeat infections typically occur in adults whose contacts are not subject to school 
term forcing. Similarly to the SIRIS model, the disease incidence/prevalence and density of 
infectives in equilibrium for the deterministic counterpart of the SIRSI model is identified 
with the class Ji. We call this scheme the repeat infection model and denote it by SIRSI. 

We finish the description of the models by pointing out that, for a fixed Po, the unforced 
{Pi = 0) deterministic counterparts of the five stochastic models introduced in this section 
have different densities of infectives in equilibrium. In Figure [2] we plot these equilibrium 
densities as a function of Pq. The values of the remaining parameters will be justified later. 
For < /3i < 0.1, the density of infectives in the deterministic versions of the five models 
oscillates with the period 1 year around the equilibrium value of the corresponding unforced 
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FIG. 2: (Color online) Infective density in the equilibrium of the unforced deterministic versions 
of the models, as a function of (Sq. The line stands for the I/N in the SIR (solid black), SIRS 
(dashed black) and SIRWS (dot-dashed black), and for the h/N in the SIRIS (solid gray) and 
SIRSI (dashed gray) model. The range of /3q = Rq{6 + fi) in the plot corresponds to i?o E [1.2, 24], 
/J, = 1/50 year~^ and 6 = 365/22 year~^. Here Rq is the so-called basic reproductive ratio of a 
disease defined as the average number of secondary cases generated by one infected individual in 
a fully susceptible population during one infectious period 



47[. The red long-dashed vertical line 



indicates the value of 00 for which Rq = 17. Parameters: 1/7 = 20 years (SIRS model); I/7 = 10 



years (SIRWS model, 



nil); 1/7 = 20 years, tj = 0.5, and a = 0.25 (SIRIS model. 



my, 1/7 = 20 



years, ij = 0.1, and a = 0.3 (SIRSI model). 

equations (see the Supplementary Material). Note that the unforced SIRWS model may 
have also limit cycles as stable attractors for some parameter values. This is, however, not 
the case for the parameter values that we will use in our analysis (throughout the main text 
/3o is fixed at the dashed vertical line in Fig. [2]). 

B. Data records 



In this paper, we consider the historical data records for pre-vaccination pertussis inci- 
dence in Greater London in the period 1946-1957 (Figure [3l the left panel) and in Ontario in 
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London records 



4.5 X 10^ 




t (year) , , , , , 

/ (llyear) 

FIG. 3: (Color online) Left: new weekly cases of pertussis in (Greater) London before vaccination. 
The dashed horizontal line is the average number of weekly cases. Right: spectrum of the time 
series. 

Ontario records {historical) 
1400 1 I 




1915 1920 1925 1930 1935 1940 
t {year) 



FIG. 4: (Color online) Left: new monthly cases of pertussis in the Canadian province Ontario before 
vaccination (upper panel) and the detrended time series (lower panel). The dashed horizontal line 
in the lower panel is the average number of monthly cases. Right: the solid black line is the 
spectrum of the full time series shown in the lower left panel, the red (black) dashed and the gray 
dashed lines are the spectra calculated from the detrended and undetrended partial time series 
shown in the same color. 
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the period 1914-1943 (Figure HI the left panels). The time series correspond to weekly data, 
in the case of London, and monthly data in the case of Ontario. Both raw and detrended 
(but not rescaled) data are shown for Ontario, where significant changes in population size 
and reporting efficiency are apparent during most of that period. In order to analyze the 
populations involved in these two data records, we shall consider only the undetrended data 
for Ontario corresponding to the last 8 years of the 29 years period. The ratio of the popu- 
lations is close to 8/3.5, according to demographic data that set the population of Greater 
London close to 8 million and that of Ontario close to 3.5 million in the period under con- 
sideration. However, the average recorded rate of new cases in a month interval, which is 
similar for the two data sets, points to a ratio of effective populations close to 1. Differences 
in surveillance coverage and/or in reporting efficiency could explain this. For London, levels 
of pertussis reporting efficiency in the range 5%-25% have been acknowledged in the liter- 



ature 



48(1 . The fact that the effective population for Ontario is approximately of the same 



size indicates that a higher reporting efficiency compensates for the smaller population. 

It is known from the case notification data from England and Wales and other locations 
that the inter-epidemic periods for pertussis in the pre-vaccination era show significant 
multiannual structure. The multiann ual p eriods lie in the range of 2-3 years and additionally 
annual outbreaks can be observed 



18 



42 



49| . The power spectra computed the data confirm 



these observations (see the right panels in FigJS] and in Fig. H]). In the spectra plots, the 
shaded region marks the range of frequencies, / G [1/3,1/2] year^^, corresponding to the 
interepidemic periods 1// G [2,3] years. Because of limited length and resolution in time, 
the spectra have a short range and a low resolution in frequency, however this affects mostly 
the annual peak. The widths of the multiannual peaks are comparable with each other. For 
Ontario, the power spectra of the full time series and of the detrended and undetrended 
partial data sets show the same dominant frequency components. In what follows, we shall 
consider only the power spectrum of the undetrended data. 
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Epidemiological meaning 


Notation 


Value 


Per capita birth/death rate 




1/50 year^"*^ 


Average lifespan 




50 years 


Rate of recovery from infection 


6 


365/22 year-i 


Average infectious period 


1/5 


22 days 


Basic reproductive ratio 


Rq 


17t 


Average contact rate 


/3o 


282.39 year^i 


Amplitude of seasonal forcing 


/3i 


[0,0.1] year-i 


Rate of loss of naturally acquired immunity 


7 


[1/40,1/10] year-i 


Average infection-derived immunity period 


1/7 


[10,40] years 


Relative infectiousness of repeat to primary (SIRSI) 
or mild to severe infections (SIRIS) 


7] 


[0,1] 


Relative susceptibility of repeat to primary (SIRSI) 
or mild to severe infections (SIRIS) 


a 


[0,1] 


Boosting coefficient 


K 


[0.2,20] 



TABLE I: The epidemiological description of the parameters and the range of their values for 
pertussis. ^ For the SIRS and SIRSI models, the whole range of possible i?o values estimated from 
the average age at first infection is explored in the Supplementary Material. 



C. Parameter values 



The values of the demographic and epidemiological parameters that are well established 



in independent data sources are kept fixed (see, for example, |l8l . 1471]). These are the 
birth/death rate /i (or the average lifespan l//i), the rate of recovery from infection 6 (or 
average infectious period 1/6) and the average contact rate Pq. The value of the latter 
corresponds to the value of the basic reproductive ratio Rn = (3q/{6 + fi) reported in the 



literature for the two data sets considered in this study 18|. If the reported value of Rq 



comes from average age at infection data, together with the assumption of SIR dynamics, 
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then it has to be corrected for models with temporary immunity. In the Supplementary 
Material we explore the whole range of possible Rq values for the SIRS and SIRSI models 
and show that the analysis of Section HTTl below is not affected by this correction. We also give 
for the SIRWS and SIRIS models the dependence of Rq on the remaining model parameters, 
computed from average age at infection data on the assumption that the naive susceptibles 
are homogeneously mixed in the whole population. 

For the remaining parameters we either use the estimates found in the literature for a 
particular model or explore an appropriate range of possible values if such estimates are 
absent. For example, the accepted range for the duration of naturally acquired immunity 



for pertussis is 7 to 20 years 50|]. For the SIRS model, we take I/7 equal 20 to 40 years. 
On one hand, the large upper limit allows for a comparison with the SIR model. On the 
other hand, as we shall see, taking I/7 smaller than 20 years would make the prediction 



of the model worse. For the SIRIS model, we take I/7 = 20 years j45| and the relative 
infectiousness and relative susceptibility of mild infections ?7, a G [0, 1]. The prediction of 
this model for other values of I/7 can still be done using the SIRS model as one of the 
limiting cases of the SIRIS model. For the SIRWS model, we take the value of the boosting 
coefficient k = 20 and I/7 = 10 years considered in Ref. jll| and in addition we study 
the behavior of the model for I/7 = 40 years and k G [0.2, 10]. For the SIRSI model, we 
set 1/7 = 20 years and the relative infectiousness and relative susceptibility of repeat to 



primary infections ?7, o" G [0, 1] [42[. We discuss in the Supplementary Material the behavior 
of this model for other values of the duration of immunity. 

For practical purposes we use only three values of the amplitude of seasonal forcing 
/3i = 0,0.05,0.1. We will however discuss the predictions of the models for intermediate 
values of Pi whenever necessary. The summary of the parameter values is given in Table HI 

D. Stochastic simulations 

To study the behavior of the five candidates for the description of pre- vaccination pertussis 
we use stochastic simulations of the processes described in Fig. [T] The simulations are 
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based on a modification of Gillespie's algorithm 5l| which accounts for the explicit time 



dependence in the contact rate 



52I ]. 53| . In Sections IIII AtlTTTEt each simulation run starts 



from a random initial condition and the prevalence of the disease is recorded with a time 
step of 0.05 year for 450 years after 50 years of transient. From each simulation run the 
power spectrum of the prevalence time series is computed with the use of the discrete Fourier 
transform. The final spectra are averages of 10"^ simulations. The population size, A^, used 
in the simulations in Sections IIII AtHTTEl is 5 x 10^. 



III. RESULTS 



We investigate the dynamics of the stochastic seasonally forced SIR, SIRS, SIRWS, SIRIS 
and SIRSI models by comparing the power spectra of long time series of the five models 
in the relevant parameter ranges. We first describe the performance of the SIR model and 
then assess whether each of its four extensions improves or worsens the results for the SIR 
model. In all cases, the spectra correspond to the fluctuations around the equilibrium of 
the unforced deterministic versions of each model = 0) or around the associated period 
one year stable attractor > 0). Numerical and simulational results for the parameter 
ranges we have explored never showed evidence of other attractors. A direct quantitative 
comparison of these power spectra with the spectra shown in Fig. [3] and Fig. H] based on 
the amplitude of the multiannual peaks and their precise location is undermined by the low 
resolution in frequency and poor statistics of the two data sets, which correspond to much 
shorter periods than the ones that can be used in the numerical simulations. That is why 
our first criterion in making the comparison between the models' predictions and the spectra 
computed from the data will be the position of the multiannual peak. In order to be able 
to explain historical data records of pertussis incidence a model's spectrum should have the 
multiannual peak located in the range of frequencies, / G [1/3, 1/2] year~^ (see Fig. [3] and 
Fig. H]). This region is the shaded region shown in all spectra plots throughout the paper. 

As pointed out before, in Sections IIII AlUTTEl we use a typical population size A^ = 5 x 10^ 
in all models and compute spectra from the prevalence time series. Changing population 
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FIG. 5: (Color online) The average power spectra from simulations of the SIR model for different 
amplitudes of seasonal forcing, f3i. Parameters: /3i = [green (light gray)], 0.05 [orange (dark 
gray)], 0.1 (black). The multiannual peaks of the spectra computed from the London and Ontario 
data sets are located in the shaded region. Henceforth, all simulation spectrum plots will be 
shown in lin-log scale and the dashed line will indicate the dominant frequency of the stochastic 
multiannual peak found from the simulations. 



size and/or computing the incidence in a given interval instead of the prevalence of the 
disease may change somewhat the amplitude of the multiannual and annual peaks but does 
not influence their position. 



A. SIR model 



The model's power spectrum as a function of frequency for different values of the forc- 
ing amplitude Pi is shown in Fig. |5l The power spectra computed from the stochastic 
simulations exhibit two types of peaks: a well defined dominant annual peak due to the 
deterministic annual limit cycle and a subdominant broad multiannual stochastic peak with 
the shape and the main frequency similar to that of the unforced case. The increase in /3i 
results in a more enhanced annual peak which means that the contribution of annual epi- 
demics in the time series increases as (3i increases. The multiannual stochastic peak reflects 
the presence of noisy oscillations in the incidence time series with that dominant frequency. 
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The mechanisms 
in several papers 



w which internal noise excites these resonant fluctuations has been treated 



27 



28 



31 



32 



43[. They can be described analytically using van Kampen's 
expansion of the master equation of the underlying stochastic process around the attractor 
of the deterministic system 5J]. We point out that in this study the dominant frequency 
of the main stochastic peak in the seasonally forced models can be computed from the un- 
forced model whose analysis is easier (see Supplementary Material for more information on 
the analytical computation of the power spectra). 

For the SIR model, the frequency and the shape of the dominant stochastic peak are 
largely independent of and the peak's frequency lies in the shaded region (see the Sup- 
plementary Material for the analytical expression of the spectrum). The dominant period 
of stochastic fluctuations in the SIR model is 2.7 years. In Sections IIII BIHTTEI we will dis- 
cuss whether the SIRS, SIRWS, SIRIS and SIRSI extensions of the SIR model improve or 
worsen its performance. In particular, we will be interested to know how the position and 
the amplitude of the stochastic peak change with respect to that of the SIR model and the 
robustness of the behavior of each model under the variation of free parameters. 



B. SIRS model 



In the SIRS model, the rate of immunity waning 7 is a free parameter. In the case 
of the lifelong immunity, 7 = 0, the SIRS model reduces to the SIR model considered 
in the previous section. Figure |6] shows the model's spectrum for different values of the 
forcing amplitude (3i and two typical values of 7. In the unforced SIRS model, /3i = 0, 
the spectrum has a pronounced stochastic peak (see the Supplementary Material for the 
analytical expression of the spectrum). However, it is situated outside the region of interest 
both for 1/7 = 40 years and for I/7 = 20 years. For fixed 7 and increasing the frequency 
of the stochastic peak is slightly shifted towards the target region but it remains outside 
it even for large values of 7 and, at the same time, the power associated with the annual 
peak becomes very large. Note that if Pi is decreased to take into account that non primary 
infections should not be subject to seasonal forcing in this model the dominant frequency 
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FIG. 6: (Color online) The average power spectra from simulations of the SIRS model for different 
amplitudes of seasonal forcing, and two values of the immunity waning rate, 7. Parameters: 
1/7 = 20 years (black), 40 years [orange (gray)]. The multiannual peaks of the spectra computed 
from the London and Ontario data sets are located in the shaded region. The dashed lines indicate 
the dominant frequencies of the stochastic multiannual peaks found from the simulations. 

of the stochastic peak is moved even further away from the shaded region. We conclude by 
noting that the SIRS model's spectrum is incompatible with the London and Ontario data 
unless in the limit of 7 — t- when it approaches the SIR model. The performance of the 
SIRS model for other values of Rq is also discussed in the Supplementary Material. 

C. SIRIS model 

Next we analyze the power spectra from simulations of the SIRIS model. This model has 
three free parameters apart from the amplitude of seasonal forcing Pi. the rate of immunity 
loss 7, and the relative infectiousness 77 and relative susceptibility a. For the limiting case 
of ?7 = and (7 = the model reduces to the SIRS model whose multiannual peaks are 
located outside the shaded region. The model's behavior is, however, well understood if 
we vary only r] and a. The results are shown in Fig. [7] for fixed /3i = 0.05 and I/7 = 20 
years. The spectrum for small values of 77 and a is reminiscent of the SIRS model. A 
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FIG. 7: (Color online) The average power spectra from simulations of the SIRIS model for f3i = 
0.05, 1/7 = 20 years and different values of the relative infectiousness, ij, and relative susceptibility, 
cr, of mild to severe infections. The spectra for other values of /3i are given in the Supplementary 
Material. Parameters: 77 = 0.1 (black), 0.5 [orange (dark gray)] and 0.75 [green (light gray)]. The 
values of rj and a were chosen from (0, 1) so as to guarantee that the equilibrium incidence in the 
unforced deterministic SIRIS model is of the same order of magnitude as in the SIR(S) models. 
The multiannual peaks of the spectra computed from the London and Ontario data sets are located 
in the shaded region. The dashed lines indicate the dominant frequencies of the multiannual peaks 
found from the simulations. 

comparison of the spectra in all panels of Fig. [7] shows that the stochastic and deterministic 
peaks are the highest for the smallest values of both 77 and a. This is because it is easier 
to generate incidence oscillations in this parameter range. The increase in the values of rj 
and a results in a flattening of the spectrum occurring through a gradual disappearance of, 
first, stochastic and, then, deterministic peaks. Moreover, the study of the dependence of 
the power spectrum on the amplitude of seasonal forcing Pi does not change this picture. 
For example, for the smallest values of 77 and a used in Fig. [7] we observe that the position 
and the shape of the stochastic peaks agree perfectly for Pi G [0,0.1] and that they are 
always outside of the region of interest (this result is given in the Supplementary Material). 
The parameter region with higher values of 77 and a is even less interesting because of 
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FIG. 8: (Color online) The average power spectra from simulations of the SIRWS model for different 
amplitudes of seasonal forcing, and two values of the immunity waning rate, 7. Parameters: 
1/7 = 10 years (black) and 40 years [orange (gray)]. The immune boosting coefficient k = 20 
in all panels. The spectra for other values of k are given in the Supplementary Material. The 
multiannual peaks of the spectra computed from the London and Ontario data sets are located in 
the shaded region. The dashed lines indicate the dominant frequencies of the multiannual peaks 
found from the simulations. 

the flattening of the spectrum. Thus for all relevant ranges of 77, a and (3i the stochastic 
multiannual peak in this model is situated outside the region of interest except when it 
approaches the limit of the SIR model. 

D. SIRWS model 

Here we show the power spectra computed from simulations of the SIRWS model for 
boosting coefficient k = 20 and different values of 7 and see Fig. [HI The results for 
K = 0.2, 1 and 10 are given in the Supplementary Material and will be discussed later in this 
section. In the absence of seasonality, /3i = 0, the spectra have a dominant multiannual peak 
situated in the target region for the whole range of I/7 G [20,40] years. When /3i increases 
a deterministic annual peak appears in the spectrum while the dominant frequency of the 
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multiannual peak stays unchanged. Similarly to the SIR model, the higher the amplitude 
of seasonal forcing Pi is the higher the annual peak is. However, unlike in that model, this 
phenomenon is accompanied by an overall complication of the structure of the spectrum for 
higher values of (3i. For example, the power spectrum for I/7 = 10 years and /3i = 0.1 
(black line in the right panel of Fig. |8]) has several secondary multiannual peaks, one of 
which is situated near but outside the shaded region. Thus, the SIRWS model predicts a 
rather broad range of frequencies corresponding to the stochastic fluctuations around the 
annual deterministic cycle. 

Finally, we discuss the results for k ^ 20 (see Supplementary Material for details). For 
K, = 0.2,1,10 and I/7 = 40 years, the stochastic multiannual peaks of all power spectra 
lie in the shaded region. The same happens for k = 10 and I/7 = 10 years but not for 
k = 0.2, 1. For K — > and 7^/2 the SIRWS model exhibits a SIRS-like behavior, thus the 
multiannual peaks of the SIRWS model are outside of the region of interest likewise in that 
model. 

E. SIRSI model 

We finish by presenting the results from simulations of the SIRSI model. Figure [9] shows 
the spectra for different values of f3i G [0, 0.1], a, 77 G (0, 1) and I/7 = 20 years. Similarly to 
the SIR model, the simulation spectra of the SIRSI model have a high stochastic multiannual 
peak in addition to the annual peak due to the deterministic limit cycle. The amplitude of 
the deterministic peaks and therefore the power associated with them increase as /3i increases 
while the stochastic peaks stay almost unchanged. The increase, however, is smaller in the 
SIRSI model than in the SIR model. This means that, notwithstanding the similarity of 
their simulation spectra, the contribution of annual epidemics is less in the former than in 
the latter, corresponding to a qualitative difference in the times series of the two models. 

For three sets of parameter values chosen in Fig. [9] the dominant stochastic multiannual 
peaks lie in the shaded region, moreover their frequencies and shapes are largely independent 
of (3i. To explore the consequences of varying r] and a for fixed 7, we considered in the interval 
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FIG. 9: (Color online) The power spectrum from simulations of the SIRSI model. Parameters: 
Pi = [green (light gray)], 0.05 [orange (dark gray)] and 0.1 (black), I/7 = 20 years. The values of 
the relative infectiousness r] and the relative susceptibility a of repeat to primary infections were 
chosen from a set of all possible values explored in the Supplementary Material. The multiannual 
peaks of the spectra computed from the London and Ontario data sets are located in the shaded 
region. The dashed lines indicate the dominant frequencies of the multiannual peaks found from 
the simulations. 

(0, 1) nine equally spaced values and constructed a grid of 81 points in the (77, a) space. As 
the position and the shape of the stochastic peak in this model is independent of (3i we 
can calculate the stochastic peak's frequency from the unforced model (see Supplementary 
Material for more details). For example, for I/7 = 20 years used in Fig. [HI the number 
of spectra for the SIRSI model with the stochastic peak within the shaded region is 62 
(for 1/7 G [10,40] years this number varies from 43 to 78). We thus conclude that the 
SIRSI model's spectrum illustrated in Fig. |9] is robust with respect to variation of all free 
parameters, and that it exhibits some degree of variability in the amplitude of the stochastic 
multiannual peak as parameters are varied within their accepted ranges. The robustness of 
these results with respect also to variations of Rq keeping the average age at first infection 
fixed at its value for the SIR model is illustrated in the Supplementary Material. 
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FIG. 10: (Color online) Incidence of the disease in a simulation of the SIRWS, SIR and SIRSI 
models. The horizontal dashed line is the average number of weekly (London) or monthly (Ontario) 
cases (compare these with the dashed lines in Fig. [3] and in Fig. Parameters: Pi = 0.05 (all 
plots); K = 20, 1/7 = 10 years (SIRWS model, the black hne in the middle plot of Fig. E]); I/7 = 20 
years, rj = 0.1, a = 0.1 (SIRSI model Ontario) and a = 0.3 [SIRSI model London, the orange (dark 
gray) line in the middle plot of Fig. [9]. Note that the level of seasonal forcing /3i does not change 
the average incidence. 

F. Comparison with real data 

The results for the power spectra of stochastic simulations show that out of the four 
extensions of the SIR model under consideration only SIRWS and SIRSI have frequency 
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values for the stochastic multiannual peak that preserve or improve the compatibihty with 
data for pre- vaccination pertussis previously obtained with the SIR model. Here we compare 
the other features of spectra of the SIR model and of the SIRWS and SIRSI extensions with 
the spectra for London and Ontario data sets shown in Figs. [3] and |H 

As a first step towards achieving this we determine the population size for each model. 
As discussed before both data sets have the similar average recorded rate of new cases 
in a month interval and thus correspond to a similar effective population. However, the 
deterministic counterparts of the models, and consequently the stochastic models too, predict 
different densities of infectives in the steady state. To resolve this, for each of the three 
models we make simulations with the same length and the same sampling time as those of 
the London or Ontario data sets and record the incidence of the disease in the sampling 
interval. The transient period in all simulations is taken to be 50 years and each simulation 
starts from a random initial condition. We calibrate the effective population size A^ for a 
given model by imposing that the incidence averaged over 10'^ different simulation runs is 
the same as the time averaged incidence in the corresponding data set, given by the dashed 
line in Fig. [3] and in Fig. |H In Fig. [10] we show typical incidence time series for /3i = 0.05 
and the effective population size A^ for each of the three models from which individual power 
spectra are computed. 

The second step is to determine for each model the level of seasonality that corresponds 
to the data set. We estimate the value of the amplitude of the seasonal forcing /3i as the 
value for which the amplitude of the annual deterministic peak in the power spectrum of 
incidence times series averaged over 10^ different simulation runs is equal to the amplitude 
of the annual peak in the spectrum computed from the data. 

Using the values of A^ and of f3i determined in this way, we produce from the stochastic 
models time series mimicking the London or Ontario data sets, and average their power 
spectra over 10^ simulation runs. The power spectra from simulation [dashed red (black) 
lines] and the spectra computed from the data sets [black solid lines] are shown in Fig. [HI 
The power spectra obtained with the SIR model have stochastic multiannual peaks with 



23 



SIRSl model (Ontario) 



1.8 X Iff 




SlRSl model (London) 



.-Kl(f 



4xl(f 



3x10' 



2xl(f 



ixia 



1.2 / l(f 



SIR model (Ontario) 



SIRWS model (Ontario) 




xlCP 



6xl(f 



4xl(f 



2xl(f 



0.2 0.4 0.6 0.8 1.0 
SIR model (London) 




0.0 





\ 

1 
1 


;\/ 


1 
I 
1 


n 

1 » 
1 ' > 


1 
I 
1 


1 > t 

■A' 


1 


,1 II 


i 


■III 
1',' > 

/'.' — >- 
/. 







0.2 0.4 0.6 
SIRWS model 



0.8 1.0 
(London) 



1 2.5 - 10' 




f (1/year) f (1/year) f (1/year) 



FIG. 11: (Color online) The solid black line is the power spectrum of the London or Ontario data 
sets as indicated in each panel. The dashed red (black) line is the average power spectrum of 
the incidence time series computed from 10'^ simulation runs with the same length and the same 
sampling time as those of the London or Ontario data sets. The dashed gray lines are four typical 
power spectra computed from an individual simulation. Parameters: /3i = 0.06 (SIRSI Ontario), 
/3i = 0.07 (SIRSI London), /3i = 0.05 (SIR Ontario), /5i = 0.06 (SIR London), /3i = 0.03 (SIRWS 
Ontario), /3i = 0.05 (SIRWS London). The remaining parameter values are as in Fig. [TUl 

amplitudes that are similar, though on average larger, than those of the data sets power 
spectra. For the SIRWS extension, the average amplitude of the stochastic peak is even 
larger, while the SIRSI extension gives a better agreement with this feature of the data with 
respect to the SIR model. 

The results of Fig. [TT] thus suggest that the SIRSI model is the one that best reproduces 
the stochastic properties of both data sets. However, as shown by the four examples plotted 



24 



in the figure for eacfi case (dasfied gray lines), tlie poor statistics of the short term samphngs 
translates into a considerable variability in the power spectra computed from each time series. 
The two historical data records considered here illustrate the principle that the incidence 
power spectra can be used to discriminate the three variants of the model, but a final 
conclusion on their performance would require the analysis of longer time series. 

IV. DISCUSSION AND CONCLUSIONS 

Different assumptions about disease and vaccine induced immunity have been proposed in 
the recent literature to model the dynamics of pertussis in the presence of mass vaccination. 
Our motivation has been to study separately the influence of disease acquired immunity only, 
in order to reduce the number of free parameters, and to consider explicitly the stochastic 
properties of the incidence time series as part of the model's output, in order to look for 
further constraints. 

We have used the basic SIR model and four variants that reflect different immune re- 
sponses. We characterized the dynamic patterns of each of these five models through the 
average power spectra of an ensemble of long term prevalence time series obtained from 
stochastic simulations. We have found that the power spectra of the five alternative models 
show quantitative and even qualitative differences. One of them (SIRS) fails to reproduce 
a stochastic peak in the frequency range found in real data for pertussis, and another one 
(SIRIS) is too stiff to produce stochastic fluctuations with the observed amplitudes. 

Assessing the performance of the other two variants, SIRWS and SIRSI, with respect 
to the SIR model requires a quantitative comparison of the simulated power spectra with 
real data, for which purpose the effective population size and the forcing amplitude of each 
data record must be determined. We illustrated this procedure by considering two publicly 
available historical data records for pre-vaccination pertussis incidence. Based on these 
relatively short time series, we found that the output of the SIRSI model agrees with the 
main aspects of the phenomenology of the disease in a broad parameter range. It reproduces 
the multiannual interepidemic periods and the amplitude of fluctuations found in the data 
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with a better agreement than SIR for the latter, while the SIRWS variant tends to produce 
stochastic fluctuations even larger than SIR. 

Considerations about the recovery proflle and its influence on the model's behavior concur 
in favoring models such as SIRSI, with a much flatter spectrum than SIRWS. Here we have 
followed the common approach of taking recovery as a constant rate stochastic process, 
instead of more realistic unimodal recovery proflles that have been proposed and studied in 



the literature 



30 



55| . For unforced models, it is known [30| that the effect of this change 



is twofold, on one hand a small displacement of the main stochastic peak towards higher 
frequencies, and on the other an enhancement of the fluctuations. We have checked (results 
not shown) that this holds too for the SIRSI and SIRWS models in the parameter region 
explored in the previous section when we change from an exponentially distributed recovery 
proflle to a unimodal gamma distributed recovery proflle with the same average. Therefore, 
under more realistic recovery proflles, the output of the SIRWS model would be further away 
from the target power spectra, as would that of the SIR model as well. 

However, the short length of the time series used in our illustrative study prevents any 
deflnite conclusion about the performance of the SIR model alternatives. Indeed, the ensem- 
ble of the power spectra obtained from simulated time series of this length exhibits enough 
variability to accommodate the two samples of real data. An extension of this analysis to a 
larger set of data could help to further elucidate the dynamics of disease acquired immunity, 
setting a solid ground for more complex models to deal with vaccine acquired immunity. 
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Supplementary Online Material Outline 

I. Deterministic counterparts of the SIR, SIRS, SIRIS, SIRWS and SIRSI models. 

II. SIRIS model for Pi ^ 0.05. 
II. SIRWS model for k ^ 20. 

IV. Analytical computation of the dominant frequency of the multiannual peak. 

V. Estimation of Rq based on the average age at primary infection. 

VI. Analytical power spectrum of the stochastic SIR and SIRS models for /3i = 0. 

I. DETERMINISTIC COUNTERPARTS OF THE SIR, SIRS, SIRIS, SIRWS AND 
SIRSI MODELS 

A. SIR model 

In the infinite population limit, — oo, the equations of the SIR model for the suscep- 
tible and infective densities are 



- = ^^{l-s)-ms^, (1) 
di 

- = P{t)si - {6 + fi)i, (2) 



where 

/3(t) =/3o(l + /3icos27rt). (3) 

The density of recovered individuals is thus given by 

r = 1 — s — i. (4) 

For /3i = 0, seasonal forcing is absent and the system has a trivial {i* = 0, s* = 1) and a 
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nontrivial fixed point. The latter can be expressed in terms of tlie basic reproductive ratio 
Rq of a disease defined as tlie average number of secondary cases generated by one infected 
individual in a fully susceptible population during one infectious period: 



M - 1) 

/3o 



s* = — 
Rq ' 



(5) 



where 



Rq 



6 + fi' 



(6) 



Stability analysis shows that the trivial equilibrium is stable if -Rq < 1 while the endemic 
equilibrium ([5]) is stable if i?o > 1- This means that the critical value of Rq given by Eq. 
(|6]) above which the disease is endemic is 1. 

B. SIRS model 

In the infinite population limit, the equations of the SIRS model are 



where P{t) and r are given by Eqs. ([3]) and (jl]), respectively. For /3i = the non trivial 
equilibrium values s* and i* of the model satisfy 



ds 

'di 

di 

di 



/3{t)si- {S + fi)i, 



/i (1 - s) - f3(t)si + jr, 



(7) 



(8) 



i* = 



{j + fi){6 + fi){R^-l) 



s* = — 

Rq ' 



(9) 



where Rq is given by Eq. ([6]). 
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Although the SIRS model is not usually formulated with explicit classes for the naive 
susceptibles and those that become susceptible because of immunity waning, and for infec- 
tives with primary or subsequent infections, they can easily be introduced. Considering the 
fraction si of naive susceptibles and the fraction S2 of susceptibles that have been infected 
before, S2 = s — si, Eq. ([7]) for /3i = becomes 



at 



dS2 



-(3oS2i - fiS2 + 7(1 - Si- S2-i) 



(10) 



yielding in equilibrium si = fi/{Poi* + fi). Similarly, from Eq. ([8]), the evolution equations 
for the fraction of primary infections ii and other infection episodes i2, 12 = i — H, are 



l3oSii - {6 + 



dii 

m 

^ = PoS2i -{5 + 11)12, 



(11) 



yielding in equilibrium = f3osl/{6 + fi) = fi/ [7(1 — s* — i*) + /i] for the fraction of 

primary infections. 



C. SIRWS model 



In the infinite population limit, the equations of the SIRWS model for the densities of 
individuals in the susceptible, infectious and waning classes are 



- = ^I{l-s)-ms^ + 2JW, (12) 



- = + (13) 

^ = 27(r — w) — k^qWi — fiw, (14) 



where r is the density of individuals in the recovered class 

r = l — s — i — w (15) 



and P{t) is given by Eq. 



D. SIRIS model 



In the infinite population limit, the equations of the SIRIS model for the susceptible and 
infective densities are 



ds 

'dt 

dii 

lit 
di2 

'dt 



/i (1 - s) - s + Po7]i2] + -fr, (16) 

s[P{t)ii + /3or]i2]-{6 + fi)ii, (17) 
(T/3or (ii + 7722) - (^ + /i)i2, (18) 



where 

r = l — s — ii — i2 (19) 



and /3(t) is given by Eq. If /3i = the model has a trivial equilibrium and an endemic 
equilibrium where the densities of infectives and susceptibles lie between and 1. The 
analytical computation of the endemic equilibrium values of il, and s* is straightforward 
but the expressions are lengthy and that is why we do not write them explicitly. 
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E. SIRSI model 



In the infinite population limit, the equations of the SIRSI model for the susceptible and 
infective densities are 



dsi 

It 

dii 

lit 

dS2 

m 

di2 
Hi 



si)- si[p{t)ii + Pom2]. (20) 

si[/3(t)zi + /3o^Z2] - + (21) 

7r - /iS2 - cr/3oS2 {ii + 77^2) , (22) 

cr/3o-52 (^1 + ^7^2) -{5 + /i)«2, (23) 



where 

r = I - si - ii - S2 - 12 (24) 



and /3(t) is given by Eq. ([3]). It can be shown that for /3i = the model has a trivial 
equilibrium and an endemic equilibrium where the densities of infectives and susceptibles 
lie between and 1. As for the previous model, we do not give its explicit form. 

A comparison of the equilibrium infective densities corresponding to the nontrivial fixed 
point as a function of /So for the five models is shown in Fig. 2 of the main text. 

II. SIRIS MODEL FOR A / 0.05 

In the main text we presented the results for the SIRIS model for fixed /3i = 0.05. 
Figure [1] shows the power spectra of the number of infective individuals from simulation 
for /?! = 0,0.05,0.1. The stochastic peaks agree perfectly for all /3i G [0,0.1] and they are 
always far away from the shaded region. There is a remarkable resemblance of these spectra 
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FIG. 1: (Color online) The dependence of the power spectrum of the 
SIRIS model on the amplitude of seasonal forcing f3i. Parameters: pi = 
[green (light gray)], 0.05 [orange (dark gray)], and 0.1 (black). The rest of the parameter 
values as for the black line in the left panel of Fig. 7 of the main text. 

with those of the SIRS model for the same value of immunity waning rate (compare them 
with the black lines in Fig. 6 of the main text). 



III. SIRWS MODEL FOR k / 20 

In the main text we presented the results for the SIRWS model for fixed boosting coeffi- 
cient K, = 20. Figure |2] shows the power spectra of the number of infective individuals from 
simulation for n = 0.2, 1, 10 and I/7 = 10 years (top panels), 40 years (bottom panels). In 
each panel there are three spectra corresponding to different values of the seasonal forcing 
amplitude. For given 7 and n the dominant stochastic peaks agree well for all Pi G [0,0.1]. 
We also observe that they are situated within the shaded region unless in the limit of the 
SIRS model (k is small and 7 ^ /.i). 
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FIG. 2: (Color online) The dependence of the power spectrum of the SIRWS model on the boosting 
coefficient k. Parameters: /3i = [green (light gray)], 0.05 [orange (dark gray)], and O.I (black). 

IV. ANALYTICAL COMPUTATION OF THE DOMINANT FREQUENCY OF 
THE MULTIANNUAL PEAK 



In this section we will describe how the dominant frequencies of the multiannual stochastic 
peaks in power spectra from simulations can be computed analytically. First, we note that 
in the three most interesting models which are the SIRWS, the SIR and the SIRSI model 
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these frequencies are independent of /3i G [0,0.1] (see Fig. 5, Fig. 8 and Fig. 9 of the main 
text). Thus, to predict the position of the stochastic peak in this parameter range it is 
sufficient to consider the unforced models with /3i = whose analysis is easier. The method 
we are going to use is usually referred to as the system size expansion first developed by van 
Kampen It has been successfully applied to several epidemiological models based on the 
S(E)IR or SIRS dynamics. We refer the reader to 0-0] for a more extensive discussion of 
the method in the context of unforced epidemiological models. The analysis of seasonally 
forced models can be found in [l, 6]. Here we will outline the most important steps of the 
calculations for the unforced SIRSI model because it has the largest number of classes. 

The analytical treatment of the stochastic SIRSI model simulated using Gillespie's al- 
gorithm starts with the construction of the master equation We describe the state of 
the system / = (/i, /2, ^3, ^4) by four numbers which are, sequentially, the number of sus- 
ceptible Si, infective Ji, susceptible 5*2, and infective I2 individuals. The total population 
size is fixed, that is why the number of recovered individuals R at any time is given by 
N — li — I2 — h — h- The transition rates T/' from the initial state I to the final state I' 
associated to the processes postulated in the diagram (e) of Fig. 1 of the main text can be 
written as follows: 

1. Immunity waning: 



Ktf:M''=l(N-h-k-h-h). (25) 

2. Infection: 5*1 Ji and 5*2 I2, where A'^ = /9o(A + vh)/N. 

K7M'''' = Mi{h + vh)/N, (26) 

7^:f:f:i'''^' = ^(3ok{i2+vk)/^ (27) 
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3. Recovery: /i,/2 — > R- 



= ^^2, (28) 

T^i^i'' = ^h- (29) 
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4. Birth and death: R A- Si and Si, Ji, 5*2, h A R. 



'1-11,12 — 1,13,14 

h,l2,h,l4 

q-h,l2,h — l,l4, 
h,l2,h,l4 



/iiV, 


(30) 


fill, 


(31) 




(32) 




(33) 


fih. 


(34) 



The master equation for the probabihty V{l,t) of having the system in state / at time t 
in the S.RSI mode, can be obtained by substituting Eqs. into the general form for 

the master equation 



dV{l,t) 
dt 



Y.[n.v{i',t)-Ti'v{i,t) 



(35) 



Given the initial and boundary conditions, the solution of Eq. ( 135|) is equivalent to the full 
stochastic simulation of the model. The exact analytical solution of the master equation for 
the SIRSI model is not feasible but we can study its limit when the size of the system, A^, 
becomes large [2]. The limit is formalized by representing the discrete variables Ik as a sum 
of two terms: 



lk{t) = Ndkit) + ViVxfc(t) 



A. 



(36) 



The first term is the deterministic macroscopic term of the order of and the second term 
are the fluctuations of the order of \fN around the macroscopic state. To preserve the 
notation introduced in the main text we denote the functions dkif) 

dk{t) = hm (37) 

W— 5>oo iV 
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as the densities of susceptible individuals Si{t) for /c = 1 and S2{t) for k = 3, and of infective 
individuals ii{t) for k = 2 and i2{t) for k = 4. 

Substituting Eq. (137|) into Eq. ( 135|) results in an equation for the probability distribution 
Il{x,t) of the fluctuations x = (xi, X2, X3, X4) around the deterministic trajectory in which 
terms of different orders in 1/\/N can be identified. After collecting and equating terms 
of the leading order one obtains the deterministic SIRSI equations with /3i = [see Eqs. 
(19)-(22)]. The same procedure for terms of the next to leading order gives rise to a linear 
Fokker-Plank equation in the variable x 



dU{x,t) _ ^ d{xjU{x,t)) Iv-o ^.N'9^n(x,t) 



dt 



k,j 



dxi 



dxkdxj 



k,j = l,...A- (38) 



In Eq. ( I38p . A(t) is the Jacobian of the deterministic SIRSI model with Pi = 



A(t) 



-g - f^i -PoSi -PoVSi \ 

g Posi - {6 + fi) /3or]Si 

-7 -a/3oS2 - 7 -<^g - {1 + fJ') -<7PoVS2 - 7 

cr/3oS2 ag (J/3oT]S2 - (5 + fi) J 



(39) 



and B(t) is the symmetric cross correlation matrix computed directly from the expansion 



Bit) 



( gsi + /i(l + si) -gsi 

-gsi gsi + {5 + 

7r + {ag + fi)s2 -o-gs2 



\ 



\ 











-ags2 



ags2 + (5 + /i)z2 / 



(40) 



where 



9 = Poih +Vi2) 



(41) 
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and 



Si-ti- S2- l2- 



(42) 



To compute the power spectrum of the fluctuations it is useful to use the Langevin equation 
[2! associated with Eq. which reads as 



dt 



Y,Akj{t)x,{t) + Lk{t), k,j = l,...A, (43) 



where Lk{t) are Gaussian noise terms with zero mean 

(L,(t))=0 (44) 



and with the correlator 

{L,it)L,it'))=B,,mt-t')- (45) 

Finally, since we are interested in the fluctuations in the stationary state we evaluate the 
matrices A(t) and B(t) at the endemic steady state (sj,i^, 52,^2) of the deterministic SIRSI 
equations (with /3i = 0). The power spectrum of the fluctuations for infectives Ji around the 
endemic equilibrium of the unforced SIRSI model is defined as Pi^{uj) = {\x2{uj)\^)- Figure 
[3] shows P/j as a function of frequency measured in the units oi f = 00/271 for the parameter 
values used in Fig. 9 of the main text. As before, the blue shaded region corresponds to 
multi-annual periods in the range 2 to 3 years where the multiannual peaks of the spectra 
computed from London and Ontario data are located. The frequency corresponding to the 
center of the shaded region, /c, is thus 5/12 year~^. The frequency of the maximum of the 
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FIG. 3: (Color online) Analytical power spectrum of the unforced SIRSI model for parameter 
values used in the right panel of Fig. 9 of the main text. The shaded region marks the frequencies 
between 1/3 year~^ and 1/2 year~^ and the dashed blue vertical line corresponds to the frequency 
of the maximum of the spectrum. 

spectrum, fpeak, is 0.409401 year~^. We used this value to draw the dashed blue vertical line 
in the right panel of Fig. 9 of the main text. 

As described in the main text we also studied the properties of the spectra for the unforced 
SIRSI model for other parameter values. For fixed 7, we considered in the interval (0, 1) 
nine equally spaced values and constructed a grid of 81 points in the (77,0") space. Figure H] 
shows the values of fpeak on this grid for I/7 = 20 years. We further calculated the number 
of spectra Ntarget for which 
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FIG. 4: (Color online) The frequency of the maximum of the analytical power spectrum of the 
unforced SIRSI model as a function of rj and a. The plot is made for I/7 = 20 years. These values 
were used to draw dashed vertical lines in Fig. [3] of this section and in Fig. 9 of the main text. 



FIG. 5: (Color online) The amplitude A of the peak of the analytical power spectrum of the 
unforced SIRSI model as a function of rj and a. The plot is made for I/7 = 20 years. 

The former gives the total number of spectra in the target region, and the latter gives the 
number of spectra with periods within 5% of the period corresponding to the center of the 
target region. The results are summarized in Table [T] for three different values of 7. 



SIRSI model i/3i = 0) 




0.2 



14 



7 (year ^) 


1/40 


1/20 


1/10 


-^target 


78 


62 


43 


No 


28 


20 


16 


A 


0.0602494 


0.0280422 


0.0159649 


A 


0.165027 


0.166171 


0.123987 



TABLE I: Values for the quantities used in the comparison and description of the spectra for the 
unforced SIRSI model in the {rj, a) space and different values of 7. 

Another quantity of interest in the description of power spectra is the amphtude of the 
peak A. Figure [5] shows the amphtude for the same 81 points with ri,a ^ (0, 1) and I/7 = 20 
years. The value of A is particularly interesting for the A'^o spectra with their peaks in center 
of the target region. In Table |T] we give the maximum and the minimum amplitude among 
these A'o spectra for different values of 7. The amplitude of the peak plotted in Fig. [3] of 
this section is 0.0538342, in the middle of the range of A when I/7 = 20 years. 

This study of the unforced SIRSI model is applicable to the case (3i 7^ 0, because we 
observe that the amplitude of the stochastic peak does not change for the seasonally forced 
SIRSI model. 

The conclusion of this section is that the SIRSI model's spectrum is robust with respect 
to variation of all free parameters. On one hand, there is a large region of parameter space 
where the dominant frequency of the stochastic peak is centered in the region of interest. 
On the other hand, these spectra have some variability in the peak's amplitude for different 
values of T], a and 7. 

V. ESTIMATION OF Rq BASED ON THE AVERAGE AGE AT PRIMARY IN- 
FECTION 

The results for all models presented in the main text are obtained for fixed basic repro- 
ductive ratio Rq = 17. Here we discuss the effect of relaxation of this assumption. 
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7 (year ^) 


1/40 


1/20 


1/10 


-^target 


53 


41 


30 




8.63 


6.06 


4.09 




16.16 


15.77 


15.38 



TABLE II: Number of the spectra in the target region, Ntarget^ and the minimum and the maximum 
values of Rq computed on the grid of 81 points in the (77, a) space. 

A. SIRSI model 

We first consider the SIRSI model as an example. Under the assumption of permanent 
immunity, Rq can be estimated from the average age at primary infection. A, and the average 
lifespan, l//i, as -Rq = 1/ (^^-4.) 3] • Thus the parameter values used in the main text (see Table 
I of the main text) correspond to the average age at primary infection ^ = 50/17 ~ 2.94 
years. On the other hand, using Eqs. fl2Ul) . A can be found from the relation ^ = 1/ [yU + A^] 
(for notation see Fig. 1 in the main text). Here = {(3{t)ii + /3o'7^2) is the rate of infection 
of naive individuals and fi is the death rate of individuals. Evaluating this expression in the 
absence of seasonality, (3i = 0, at the endemic fixed point gives A = 1/ [/i + /3o («i + ^2)]- 
From this equation we can obtain the average contact rate, Pq, corresponding to the average 
age at infection of naive individuals estimated as 50/17 years. To achieve this, we use 
the analytical expressions for and ^2 where we substitute the numerical values for all 
parameters but /3o and solve 50/17 = l/[/i + /3o (^1 + ^7^2)] ^^i /So- 
Using for a fixed value of 7 the same grid of 81 points in the (77, a) space we finally 
recompute the maximum and minimum values for Rq = Po/{l + f^) on this grid and the 
number of target spectra with the peaks situated in the shaded region, Ntarget- The results 
are summarized in Table II. 

From comparison of the results shown in Table I and Table II we observe that the number 
of spectra in the target region reduces slightly if the assumption of fixed Rq is relaxed. 
However, we see that even for the largest value of 7 we explored at least 37% of the spectra 
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correspond to the interepidemic periods compatible with the London and Ontario data. We 
conclude that the results presented for the model in the main text are very robust with 
respect to the parametrization of the model. 

B. SIRS model 

A similar analysis can be carried out for the SIRS model, solving for /3o the equation 
50/17 = 1/ [yU + /3o^*] for each parameter choice. For the range of average immunity period 
explored in the main text, I/7 G [10,40] years, this approach yields Rq = /3o/(7 + Z^) 
in the range [3.68,8.12] and all the stochastic peaks are shifted to the left with respect 
to those shown in Fig. 6 of the main text. More specifically, we observe that for the 
whole range of Rq the peak's frequencies in the SIRS model are equal to that in the SIR 
model, fpeak = 0.37 year~^. We conclude that under the assumption that the totally naive 
population is well mixed, the SIRS model fails to produce spectra with the stochastic peak's 
frequencies significantly different from that in the SIR model. 

C. SIRWS and SIRIS models 

We finally give the expressions for the SIRWS and the SIRIS models from where Ro 
keeping the average age at first infection fixed at its value for the SIR model can be computed 
for each parameter choice. For the former, /3o is found from 50/17 = l/[/i + /3o^*]) where 
i* is the density of infectives in equilibrium of Eqs. (|T2|) -( 1T4|) . For the latter, we solve 
50/17 = 1/ [/i + /3o (^1 + (^0, where il and ^^e the densities of individuals with 

severe and mild infections in equilibrium of Eqs. ( [T6|) -( fT8|) . 
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VI. ANALYTICAL POWER SPECTRUM OF THE STOCHASTIC SIR AND SIRS 
MODELS FOR /3i = 

The analytical power spectrum of the fluctuations for infectives / around the endemic 
equilibrium of the unforced = 0) SIR, Eqs. ([I])-®, and SIRS, Eqs. ©-([H]), models can 
be computed following the approach described in Section HVl of the Supplementary Material. 
For each model, this analytical expression can be used to predict the dominant frequency of 
the stochastic peak in the corresponding stochastic seasonally-forced model if the position 
of the peak does not depend or depends only very slightly on This is typically the 
case when /3i is not very large so that the stochastic model exhibits fluctuations around the 
period one stable attractor. 

A. SIR model 

The power spectrum for the fluctuations of infectives in the unforced SIR model is given 

by 

where Z = (D - u^f + T^u^, D = n{(3o - 5 - /i) and T = -(3ofi/{5 + ^u). Here D and T are 
the determinant and the trace of the linear approximation of Eqs. ([I])-® at the endemic 
equilibrium, Eq. ([S]). We refer the reader to Ref. js] for details on the computation of the 
SIR power spectrum. 

B. SIRS model 

The power spectrum for the fluctuations of infectives in the unforced SIRS model is given 

by 
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2{6 + fi){j + fi){f3o-6-fi) 




) 



(49) 



where Z = {D - u^f + T^w^ D = (7 + /i)(/3o - 5 - /i), T = -(7 + /i)(/3o + + 1 + ^^) 



and K = (/3o + 7)^(7 + /^) — 7(/3o — 5 — /i)(7 + 5 + /i). Here D and T are the determinant and 
the trace of the hnear approximation of Eqs. dZ])-® at the endemic equihbrium, Eq. (|9]). 
We refer the reader to Refs. Q and for details on the computation of the SIRS power 
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